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AN IMPLICIT SEMIANALYTIC NUMERICAL METHOD FOR THE 
SOLUTION OF NONEQUILIBRIUM CHEMISTRY PROBLEMS 

By R. A. Graves, Jr., P. A. Gnoffo, and R. E. Boughner 

INTRODUCTION 

Many physical phenomena are modeled by systems of linear and/or nonlinear 
ordinary differential equations (see for example references 1 to 5) which are 
defined as stiff systems when a large spread in negative eigenvalues exists. 

Such stiff systems commonly arise in nonequilibrium chemistry problems involving 
kinetic and photochemical reactions. The governing equations for these stiff 
systems are difficult to solve numerically using classical techniques because 
the error growth is rapid, and unless the equations are integrated using a very 
small time step, the results can be meaningless. To alleviate the problems 
involved with stiff systems, a great deal of effort has been expended in 
developing numerical solution techniques, both explicit and implicit, for 
stiff ordinary differential equations. References 6 to 9 review some of the 
more popular numerical methods and present the results of numerical comparisons 
between the methods. A generalized conclusion resulting from the studies of 
references 6 to 8 is that the implicit methods are more desirable because of 
their increased stability and, in some instances, significantly fewer mathe- 
matical operations. In these and other studies, a rather simple (yet 
fundamental) implicit technique was not investigated because these 
studies used a generalized equation which did not take advantage of 
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the form of the governing conservation equation for chemical species. The 
governing conservation equations for systems of chemical reactions can generally 
be written in the form of first-order ordinary differential equations. These 
equations can be solved by a simple implicit semianalytic technique which is 
derived from a quadrature solution of the governing equations. This method 
is mathematically simpler than most implicit methods and has the exponen- 
tial nature of the problem embedded in the solution. 

The objective of this paper is to present the development of the semi- 
analytic technique (SAT) and to compare its efficiency to that of several of 
the more popular methods available. 

SYMBOLS 


a,b 

general coefficients, see equation (1) 

C 1 ,C 2 

curve fit coefficients, see equation 

(5) 

C 1’ C 2 

curve fit coefficients, see equation 

(6) 

HS 

Hermit e-Simpson 


*n 

t Y n " Y (X n )]/Y(X n ) 


h 

step size 


RK4 

fourth-order Runge-Kutta 


DEQ 

Adams' fourth-order P-C 


TM 

Treanor's method 


DIFSYS 

modified midpoint rule 


TR 

Trapezoidal rule 


TR-EX 

Trapezoidal rule with extrapolation 


CAL 

Calahan's method 


LW1 

Liniger-Wil loughby - Method 1 




LW3 

SAT 

A 

Y 

n 

Y(X n ) 

5 


Liniger-Willoughby - Method 3 
semianalytic technique 
eigenvalue 
calculated value 
"exact value" 

transformed coordinate, see equation 3a 


MATHEMATICAL DEVELOPMENT 


The governing equation for the conservation of chemical species in 
nonequilibrium chemically reacting systems can generally be written in the 
form of a first-order ordinary differential equation (see ref. 1-0) : 


dY. 

* a(t)Y. = b(t) (1) 


where aft) and b(t) generally represent the loss and production rates of 
species i, respectively. The solution of equation 1 in terms of quadratures 
is: (This procedure is similar to that used in ref. 11.) 


t J+l t J+l 

-J a(t)dt J+l -/ a(t ' ) dt ' 

J + 1 J +• J . ^ ■(- 

Y. = Y. e + / b(t)e- dt 

11 [j ' 


( 2 ) 


This equation can be further simplified by introducing the transformation 

t J +i 

5 = a(t')dt' 

dC = -a(t)dt 


hence, equation 2 becomes: 



■ Y i *' 51 * t 1 1§ 

where 

/** 

{, = j a(t)dt 
t J 

For the least complicated case, the coefficients a(t) and b(£)/a(0 
can be approximated by linear functions. 


a(t) = C 1 + C 2 (t - t J ) 


C : = a(t J ) 


„ _ a(t J+1 ) - a(t J ) 

2 At 


C 1 + C 2^ 


1 ' a(0) 


b(€i) b CO) 


r - aUl)a(O) 


It should be noted that due to the transformation from t to 5 that 


b(0) _ b(t J+1 ) bfCil bft J l 

aCO)-^ -d — - ^ 


Introducing equation 6 into equation 3 results in: 


f f 1 - ^ + Z 1 ( c i * C 2 «e‘ 5 d5 
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This equation can now be integrated by parts to obtain the following semi- 
analytic implicit result: 

Y'J = e + CC X + C 2 )(l - e *) - C^e - (8) 

(Note: To have stability and accuracy, it is necessary that e < 1.) 

Equation 8 is semianalytic in nature and includes the inherent exponential 
behavior of the stiff problem directly in the solution. Equation 8 must be 
solved implicitly (iteratively) as the constants £ , and C depend on the 
conditions at the advanced time t^ + ^. 

An error analysis was performed, using the method of chapter 2 of refer- 
ence 12, to determine the errors incurred in making the linear approximations 
for the coefficients a(t) and b(£)/a(£)* The lowest order error terms are: 

E = - {AeV(u) + AtVtnHCj + Y?)J. 

where y"(w) is the second derivative of the ratio b(£)/a(£) on the interval’ 
0 < w < and (5 n (n) is the second derivative of a(t) on the interval 

t J < n < t J+1 . 

Numerical Experiments : 

System I (ref. 8) 

Y' = -0.04Y 1 + 10 4 Y 2 Y 3 

Y 2 = 0,04Y 1 - 10 4 Y 2 Y 3 - 3 x 10 7 Y 2 

Y 3 = 3 x 10 7 Y 2 

Y^O) = 1 Y 2 (0) = 0 Y 3 (0) = 0 

This system is nonlinear, and no exact solution for this system exists. 

The eigenvalues, determined from the Jacobian matrix of the system at X = 0 
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are ^2 a X 2 - 0 and X^ = -0.04. * max changes from 0.04 to 2405 for 
0 < X < 0.02. The eigenvalues for 0 < X < 40 are all strictly negative or 
zero, with X^ = 0, X 2 - -10 and X^ ~ -10^ to -10 4 . The sharp increase in 

the magnitude of X, makes this a particularly difficult stiff system to 
work with. In addition this system presents some starting problems for SAT 
since a 2 (0) = 0 and hence b 2 (5 1 )/a 2 (C 1 ) is meaningless. To circumvent this 
problem for Y 2 , two techniques using constant h were tried: first, the 
Hermit e-Simp son method, reference 13, and secondly, the Runge-Kutta fourth-order 
method. The RK4 start gave the best results. Table I gives the results for 
this system on the CDC 6600 as well as the results of Lapidus and Seinfeld, 
reference 8, for the IBM 7094. (The CDC 6600 is approximately 10 times faster 
than the IBM 7094.) It should be noted that due to the nature of this system 
Y 3 (X) was calculated by Yj(X) = 1 - Y 2 (X) - Yj. (X) . 

The most successful application of the SAT was to use the RK4 one step 
to obtain Y^ (0.0005), Y 2 (0.0005), and Y^ (0.0005) and then use the SAT with 

I 

a step size of 0.2 to compute the solution from 5 X 10" 4 < X < 40. As can be 
seen in figures 1 and 2 the linear approximation for b(£)/a(£) is very 
accurate for this system. 

System II (ref. 8) 

Y* = -200 (Y - F (X) ) + F' (X) 

Y (0) = 10 

F(X) = 10 - (10 + X)e _X 
Exact Solution Y(X) = F(X) + 10e" 200X 
Here, F(X) is a slowly decaying solution component and 10e~ 20 ^ X decays 
rapidly. The large negative eigenvalue of -200 makes exp(-200X) negligible 
compared to exp (-X) in the F(X) component. Results using SAT on the 
CDC 6600 are compared to results obtained by Lapidus and Seinfeld in Table 2. 
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At equivalent step sizes, SAT produced less than or equal to the error 

encountered using other methods, and worked faster than any other method (times 
for this system are based on computing over a range 0 < X < 15) 

System III (ref. 8) 

YJ = -O.lYj - 49.9Y 2 

Y* = -S0Y 2 

Y^ = 70Y 2 - 120Y 3 

Y^O) = 2 Y 2 (0) = 1 Y 3 (0) = 2 

, Exact Solution: 


YjCX) = e 
Y 2 (X) = e 
Y 3 (X) = e 


-0.1X 


-sox 

-sox 


+ 


+ 


-50X 

e 


-120X 

e 


Eigenvalues X^ = -120,, = -50, ~ -0.1 

Because the solution components due to X^ and X 2 decay rapidly, a 
stiff method which was not restricted by the magnitude of these values is 
desired. The SAT (which yields the exact solution of Y 2 (X) for any h since 
the linear approximation to a(X) and b/a(£) gives the true variation 
of these functions) is compared to results obtained by Lapidus and Seinfeld in 
Table 3. An h = 0.01 produced results which were better than the results 
obtained by any other method. However as the step size increased, the 
accuracy dropped off rapidly and at an h = 0.2, the solution was very different 
from the exact solution (except, of course, for Y 2 (X) which remained exact). 
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After examining the problem, it was found that on an interval 
of 0 - X - 0.2, with h = 0.2, the linear approximation to b/a(£) was 
very poor. For example, b/a(X) = -400e and b^/a^(X) = (70/120) e 

(see fig. 3). The effect of these terms on b/a(£) decays rapidly after 
X = 0.2, and they can be approximated by a linear function, so a method was 
tried using h = 0.01 to arrive at X = 0.2 and then proceed from X = 0.2 
with h = 0.2. This method was the fastest and yielded reasonable results. 

System IV 

Y 1 = 0,8Y 2 ‘ °' 01Y i ‘ 1()7 y i Y 2 Y 3 + 10Y l Y 3 " 100 Y i Y 2 

= -0.8Y 2 - lOY^ + 10 6 Y 2 Y 4 +10 4 Y Y 4 
Yj = O.OIY^ + 10 7 Y 1 Y 2 Y 3 + 2000Y 4 - 10^^ 

Y 4 = -10 6 Y 2 Y 4 + lOOY^ - 20000Y 4 

Yj (0) = 0.9; Y 2 (0) = 0.05; Yj(0) = 0.05; Y 4 (0) = 0 

No exact solution for this nonlinear system was obtained* The eigenvalues 
for this system, calculated from the Jacobian using values of Y from RK4 , 
are widely separated in magnitudes* All of the eigenvalues are negative or 
0 on a range 6*146 x 10 < X < 7.36 x 10 Typical values on the range are 

at X = 7 x 10 -6 , Aj = -1.017 x 10 S , A 2 = -4.979 x 10 4 , A 3 = -2.7102 x 10 1 and 

_9 / 

^4 ~ “2.515 x 10 . Results for this system appear in Table 4. This system 

was only stiff for a short time and none of the methods had problems with 
stability on a range 0 < X < 2 x 10“5. Using RK4 with h = 1 X 10 -8 as a 
standard of comparison, the tables indicate that for any given step size, RK4 
was more accurate than any implicit method. SAT gave accuracy comparable to 
other implicit methods and ran at approximately the same speeds as these 
methods for equivalent step sizes. In this system, SAT showed no advantage 
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over any other method. However, this system does indicate that SAT gives 
reasonable results for nonstiff nonlinear systems. 

CONCLUDING REMARKS 

As developed herein the crucial approximation is the linearization of 
the coefficients within the integral in the quadrature form which allows the 
semianalytic form to be obtained. In some cases this approximation is very 
good but in some applications, the linear approximation can be in error for 
what appear to be not unreasonable time steps. As demonstrated, this problem 
was overcome by having a variable time step which is small in the region where 

the Linear Approximation is in error. 

Additionally, quadratic and exponential curve fits were tried. However 
these approximations produced results which were approximately equivalent to 
those obtained with the linear approximation. Because the linear approxima- 
tion was the simplest to program and because of the consistently good results 
it yielded, it was chosen as the most desirable approximation evaluated. 

An important feature of the semianalytic technique is that it will allow 
the computation of nonequilibrium chemical systems to and including the equi- 
librium state. For systems where the rates. are large (typical of approaching 
equilibrium) the SAT equilibrium condition is the exact solution for 
equilibrium. 

As demonstrated in the example problems, the semianalytic technique is 
both rapid and accurate and should be applicable to those stiff problems which 
can be modeled by an equation like that used in this development. 
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TABLE I.- COMPARISON OF RESULTS FOR SYSTEM I 



t. 

R] 

'N 

R 2 N 

R 3 

’N 

Time, 

sec 

Method 

h 

X = 0.4 

X = 10 

X = 0.4 

X = 10 

X = 0.4 

X = 10 

IBM 7094 

CDC 6600 

RK4 

0.001 

0 

0 

0 


0 


- 

* 

DEQ 

0.001 

- 

- 

- 

- 

- 

- 

- 


TM 

0.01 

- 

- 

- 

- 

- 

- 

- 

- 

DIFSYS 

0.001 

- 

- 


- 

- 

- 

- 


TR 

0.2 

1.35x10-3 

1.05x10-3 

2.12x10-! 

2.4x10-! 

9x10" 2 

1.5xl0- 2 

9.3 

- 

TR-EX 

0.2 

1.72x10-5 

3.6x10-4 

3.5x10-2 

4.3x10-4 

6.8x10-4 

1.2xl0- 3 

34 

- 

CAL 

0.005/0.02 

2/4x10-3 

1.01x10-1 

2.5x10° 

6.0x10"! 

1.62x10-1 

5.4x10-1 

10 

- 

LWI 

0.2 

1.6xl0- 4 

4.9x40- 4 

2.4xl0- 4 

1.3x10-4 

3.2x10-3 

4.4x10-4 

20 

- 

LW3 

0.2 

5 . 9xl0 -4 

7.1xl0- 5 

2.9x10-3 

1.1x10-3 

4xl0- 2 

1.9x10-3 

23.3 


RK4 

0.0005 

0 

O 

0 

0 

0 

0 

- 

20.7 

HS 

0.001 

2.5xl0- 5 

6.88x10-3 

1.26x10-4 

2.86xl0- 2 

1.66xl0 -3 

3.65xl0- 2 

- 

168 

SAT 

(HS Start) 

0.001 

4.07xl0- 5 

2.54xl0- 4 

1.25x10-3 

1.97x10-4 

1.87xl0- 2 

6.52xl0- 4 


37 

SAT 

(HS Start) 

0.001/0.1 

3.79xl0- 5 

2.42xl0- 4 

1.25x10-3 

1.95x10-4 

1.88xl0- 2 

6.33x10-4 

- 

4 

SAT 

(RK4 Start) 

0.0005/0.2 

2.26xl0 -6 

1.31x10-5 

1 . 62xl0 -7 

1.74x10-3 

1.97xl0- 5 

4.44xl0” 5 

- 

1.43 

SAT 

(RK4 Start) 

0.0005 to X* 0.1 

then 0.3 

' 

1.97xl0~ 5 

3.96xl0“ 5 

9.04xl0’ 5 

1.62x10-4 

1.31xl0" 3 

2.1xl0- 4 

- 

.9 








TABLE 2.- COMPARISON OF RESULTS FOR SYSTEM II 


Method 

h 

X = 0.4 

X = 10 

IBM 7094 
Time, sec 

CDC 6600 
Time, sec 

RK4 

0.01 

1.0 x 10" 5 

2.0 x 10" 9 

11 

- 

DEQ 

0.005 

3.0 x 10' 9 

2.0 x 10 _9 

18 

- 

TM 

0.2 

6.7 x 10” 8 

1.0 x 10" 9 

16.5 

- 

DIFSYS 

0.1 

(a) 

(a) 

(a) 

- 

TR 

0.2 

1.85 x 10" 2 

4.3 x 10~ 5 

2 

- 

TR-EX 

0.2 

1.4 x 10~ 4 

1.0 x 10' 8 

36 

- 

CAL 

0.01/0.2 

1.7 x 10“ 2 

4.0 x 10" 8 

1 

- 

LW1 

0.2 

1.1 x 10" 3 

5.0 x 10" 8 

3 

- 

LW3 

0.2 

1.8 x 10‘ 3 

9.0 x 10~ 8 

4 

- 

SAT 

0.2 

9.35 x 10" 4 

4.1 x 10" 8 

- 

0.09 


(a) Unstable 








TABLE 3.- COMPARISON OF RESULTS FOR SYSTEM III 


Method 

h 

R 


R 2 

Z N 

R 3 

•*N 

Time, sec 

Time, sec 


X = 0.4 

X = 10 

X = 0.4 

X = 0.4 

IBM 7094 


RK4 

0.01 

2.0 x 10-7 

5.4 x 10-7 

3.0 x 10" 1 

3.0 x 10 _1 



deq 

0.01 

2.0 x 10" 4 

8.1 x 10" 7 

9.5 x 10 _1 

7.4 x 10 5 

23 : 

- 

TM 

0.2 

4.0 x 10 -4 

1.35 x 10~ 4 

1.1 x 10 5 

1.2 x 10 5 

1 

- 

DIFSYS 

0.1 

5.0 x 10‘ 4 

2.16 x 10 -4 

9.4 x 10 1 

8.3 x 10 2 

22 

- 

TR 

0.2 

1.0 x 10‘ 3 

-2.7 x 10" 4 

6.5 x 10 7 

1.3 x 10 5 

1.3 

- 

TR-EX • 

0.2 

4.0 x 10 -5 

8.1 x 10" 4 

5.7 x 10 1 

8.0 x 10 1 

30 

- 

CAL 

0.01/0.2 

2.0 x 10' 3 

2.7 x 10" 6 

2.5 x 10 5 

1.6 x 10 5 

1 

- 

LW1 

SAT 

0.2 

4.0 x 10~ 3 

1.1 x 10" 2 

5.0 x 10 5 

5.0 x 10 5 

3 

- 

0.01 to X — 0.2 
0.2 after X = 0. 2 


2.13 x 10~ 2 

2.08 x 10" 2 

0 

5.3 x 10 2 

- 

0.27 












TABLE 4.- COMPARISON OF RESULTS FOR SYSTEM IV 


Method 

— 

h 

Time, sec 

R 


R 

2 N 

R 

5 n 

R 4 

N 

CDC 6600 

QQRQj 

X=2xl0~ 5 

X=lxl0" 6 

X=lxl0~ 6 

X=lxl0 -6 

X=2xl0“ 5 

X=lxlO -6 

X=2xlO -5 

RK4 

10‘ 8 

2.94 

0 

0 

0 

0 

0 

0 

0 

0 

RK4 

10" 6 

0.0584 

9.17xl0‘ 6 

3.11xl0~ 3 

0 

0 

9.24xl0" 3 

4.21xl0~ 6 

1.38xl0 -5 

1.26xl0" 4 

TRAP 

2xl(T 8 

6.11 

0 

0 

0 

8xl0~ 6 

2.56xl0‘ 6 

0 

2.31xl0" 6 

0 

TRAP 

icf 6 

0.177 

4.67xl0' 4 

7. lxl O -2 

0 

l.OxlO" 5 

5. 24xl0 -3 

9.59X10' 5 

9. 73x10 4 

8.9xl0" 5 

SAT 

5x1 0 -8 

4.3 

1.94xl0' 5 

5.47xl0" 3 

5.90xl0" 4 

5. 82xl0~ 4 

2.04x10 4 

9.48xl0' 6 

4.48xl0' 4 

5 . 71xl0~ 4 

SAT 

10" 6 

0.159 

4. 65x1 0 -4 

9.62xl0 -3 

1.6xl0‘ 5 

6. OxlO -6 

5.75xl0' 4 

8.98xl0" 4 

1.30xl0 -3 

1.32xl0 -3 































